#include <istream>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <pthread.h>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <stdlib.h>

#include "UPGMA_nearest_neighbor.h"
#include "KITSCH_nearest_neighbor.h"
#include "bcpp_tree.h"
#include "bcpp_tree_ops.h"
#include "portable_timer.h"
#include "perlutil.h"
#include "matrix.h"
#include "get_arg.h"
#define DEBUG_LEVEL 2
#include "debug_func.h"

using boost::format;
using boost::lexical_cast;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::cerr;
using std::ifstream;
using std::map;
using std::pair;
using std::istream;
using std::max;
using std::min;


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	print and analyse orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
bool open_filename(ofstream& ofs, const string& file_name)
{
	ofs.open(file_name.c_str());
	if (ofs)
		return true;
	cerr << "Error\nCould not open " << file_name << "\n";
	return false;
}
bool print_gene_set(ofstream& ofs, const string& file_prefix, const string& file_name,
							const string& file_suffix,
							const string& header, const bcpp_tree_ops::t_ortholog_data& data)
{
	if (!data.gene_sets.size())
		return true;

	if (!ofs.is_open() &&
		!open_filename(ofs, file_prefix + file_name + file_suffix))
		return false;

	cout << boost::format("%|2$+3| sets from %|1$-10|:\n"
							"\t%|3$+5| tax1 genes\n"
							"\t%|4$+5| tax2 genes\n")
				% file_name
				% data.gene_sets.size()
				% data.cnt_tax1
				% data.cnt_tax2;
	cout.flush();
	ofs << header << std::endl;
	for (unsigned i = 0; i < data.gene_sets.size(); ++i)
	{
		for (unsigned j = 0 ; j < data.gene_sets[i].size(); ++j)
			ofs << data.gene_sets[i][j].name << ":"
				<< data.gene_sets[i][j].len << "\t";
		ofs << "\n";
	}
	ofs << std::endl;
	return true;
}


struct t_outputstreams
{

	t_outputstreams(const string& prefix) : file_prefix(prefix)
	{
		open();

		ofs_all_nodes	<< "To leaf\tTo root\tGene id\tNode type\n"
						<< string(80, '-') << "\n";

		ofs_full_tree	<< "[Branch length / To root / To leaf] Gene id {Node type}\n"
						<< string(80, '-') << "\n";
		ofs_timing <<	"Genes\tSum of squares\tAverage %stddev\tCount Trees\tbest iteration\tTime\n"
						"-----\t--------------\t---------------\t-----------\t--------------\t---------\n";
	}
	string	 file_prefix;
	ofstream ofs_timing;
	ofstream ofs_ortholog_sets[bcpp_tree_ops::do_collect_orthologs::eORTH_TYPE_CNT];
	ofstream ofs_full_tree;
	ofstream ofs_newhampshire_tree;
	ofstream ofs_all_nodes;
	bool open();
	bool printout_analysis(	bcpp_tree::t_tree&		tree,
							const string&				header,
							bcpp_tree_ops::do_calculate_branch_lengths&
														tree_brach_lengths,
							const string&				tax1,
							const string&				tax2);
};


bool t_outputstreams::printout_analysis(bcpp_tree::t_tree& tree, const string& header,
										bcpp_tree_ops::do_calculate_branch_lengths&
															tree_branch_lengths,
										const string& tax1, const string& tax2)
{
	// set type for each node
	bcpp_tree_ops::do_set_node_types tree_orthology(tree, tax1, tax2);



	//	cerr << "\n" << string(80, '*') << "\n\n";

	// dump all nodes
	{
		ofs_all_nodes << header << std::endl;
		bcpp_tree_ops::op_print_tree dump_all_nodes
											(tree,
											ofs_all_nodes,
											tree_branch_lengths.lengths_to_leaves,
											tree_branch_lengths.lengths_to_root,
											tree_orthology.node_types,
											true // nodes only
											);
		ofs_all_nodes << std::endl;
	}

	// print full tree
	{

		ofs_full_tree << header << std::endl;
		bcpp_tree_ops::op_print_tree print_tree
											(tree,
											ofs_full_tree,
											tree_branch_lengths.lengths_to_leaves,
											tree_branch_lengths.lengths_to_root,
											tree_orthology.node_types,
											false // nodes only
											);
		ofs_full_tree << std::endl;

		ofs_full_tree << header << std::endl;
		bcpp_tree_ops::op_print_tree_elaborate print_tree_elaborate
											(tree,
											ofs_full_tree,
											tree_branch_lengths.lengths_to_leaves,
											tree_branch_lengths.lengths_to_root,
											tree_orthology.node_types,
											3	// precision
											);
		ofs_full_tree << std::endl;
	}

	// print new_hampshire tree
	{
		ofs_newhampshire_tree << header << std::endl;
		bcpp_tree_ops::do_newhampshire_format_str nh_str(ofs_newhampshire_tree);
		tree.for_each(nh_str);
		ofs_newhampshire_tree << std::endl;
	}

	// print gene sets
	{
		bcpp_tree_ops::do_collect_orthologs data(tree_orthology.node_types);
		tree.for_each(data);
		std::string file_names[bcpp_tree_ops::do_collect_orthologs::eORTH_TYPE_CNT] =
											{
													"ortho_1_to_1"    ,
													"ortho_1_to_m"    ,
													"ortho_m_to_1"    ,
													"ortho_m_to_m"    ,
													"orphaned_1_to_1" ,
													"orphaned_1_to_m" ,
													"orphaned_m_to_1" ,
													"orphaned_m_to_m" ,
													"orphaned_tax1"   ,
													"orphaned_tax2"  };
		for (unsigned i = 0; i < bcpp_tree_ops::do_collect_orthologs::eORTH_TYPE_CNT; ++i)
		{
			if (!print_gene_set(ofs_ortholog_sets[i], file_prefix, file_names[i], ".data",
													header, data.ortholog_sets[i]))
				return false;
		}
	}
	return true;

}



bool t_outputstreams::open()
{
	return(
			open_filename(ofs_full_tree, 		file_prefix + "full_topology.tree"	)	&&
			open_filename(ofs_newhampshire_tree,file_prefix + "new_hampshire.tree"	)	&&
			open_filename(ofs_all_nodes, 		file_prefix + "all_nodes.list"		)	&&
			open_filename(ofs_timing, 			file_prefix + "timing"		)

		);
}

struct t_incorrect_pairs
{
	string cluster_name;
	unsigned cnt_genes;
	unsigned cnt_expected_pairs;
	unsigned cnt_missing_pairs;
	unsigned cnt_duplicate_pairs;
	t_incorrect_pairs(const t_matrix_io& data, const string& cluster_name_):
		cluster_name(cluster_name_),
		cnt_genes(data.cnt_names),
		cnt_missing_pairs(data.cnt_missing_pairs),
	cnt_duplicate_pairs(data.cnt_duplicate_pairs)
	{cnt_expected_pairs = cnt_genes * (cnt_genes - 1) / 2;}
};


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main (int argc, char *argv[])
{
	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;


	//
	// open input file
	//
	istream* datafile = &cin;
	ifstream infile;
	if (args.ks_file.length())
	{
		infile.open(args.ks_file.c_str());
		if (!infile)
		{
			std::cerr << "Error opening [" << args.ks_file << "] \n";
			return 1;
		}
		datafile = &infile;
	}

	//
	//	prepare files for output
	//
	t_outputstreams strms (args.out_file_prefix);


	unsigned cnt_cluster = 0;
	//
	//	save clusters with missing entries
	//
	vector<t_incorrect_pairs> problems;



	{
		while (datafile->good())
		{
			++cnt_cluster;


			//
			//	read data from input file
			//
			vector<t_matrix_io>	data;
			string				cluster_name;

			args.VERBOSE() << "\tReading gene pairs...\n";
			{
				t_progress_indicator dots(args.VERBOSE(), 20000);
				read_matrix_from_stream_in_fasta_pairs_format(*datafile, dots, data, cluster_name);
			}



			if (!data.size())
				break;

			if (!cluster_name.length())
				cluster_name = lexical_cast<string>(cnt_cluster);
			for (unsigned i = 0; i < data.size(); ++i)
			{
				string sub_cluster_name = cluster_name;
				if (data.size() >= 0)
				{
					sub_cluster_name += '.' + lexical_cast<string>(i);
				}

				t_matrix_io& sub_cluster = data[i];
				// save problems
				if (sub_cluster.cnt_duplicate_pairs || sub_cluster.cnt_missing_pairs)
					problems.push_back(t_incorrect_pairs(sub_cluster, sub_cluster_name));

				//
				//	write out number of genes
				//
				strms.ofs_timing << boost::format
											("%|1$-5|\t")
											% sub_cluster.cnt_names;
				strms.ofs_timing.flush();


				double sum_sqrs				= 0.0;
				double avg_percent_stddev	= 0.0;
				unsigned cnt_trees_examined = 1;
				unsigned best_iteration		= 0;
				t_portable_timer				timer;
				bcpp_tree::t_tree paralog_tree;

				//888888888888888888888888888888888888888888888888888888888888888888888888

				// UPGMA

				//888888888888888888888888888888888888888888888888888888888888888888888888
				if (args.analysis_type == eUPGMA)
				{
					args.VERBOSE() << "\tBuilding UPGMA tree...\n";
					t_progress_indicator dots(args.VERBOSE(), 50000);
					upgma_tree::upgma_make_tree(paralog_tree, sub_cluster.names, sub_cluster.matrix,
												args.fill_lacunae, dots);

					// calculate branch lengths
					bcpp_tree_ops::do_calculate_branch_lengths branch_lengths(paralog_tree);

				}
				//888888888888888888888888888888888888888888888888888888888888888888888888

				// kitsch optimized UPGMA

				//888888888888888888888888888888888888888888888888888888888888888888888888
				else if (args.analysis_type == eUPGMA_KITSCH)
				{
					args.VERBOSE() << "\tBuilding initial UPGMA tree...\n";
					{
						t_progress_indicator dots(args.VERBOSE(), 50);
						upgma_tree::upgma_make_tree(paralog_tree, sub_cluster.names,
													sub_cluster.matrix,
													args.fill_lacunae, dots);
					}

					// calculate branch lengths
					bcpp_tree_ops::do_calculate_branch_lengths branch_lengths(paralog_tree);

					// get sum squares for upgma part
					bcpp_tree_ops::do_sum_sqrs2 analyse_tree(sub_cluster.matrix,
//					bcpp_tree_ops::do_sum_sqrs analyse_tree(sub_cluster.matrix,
																branch_lengths.lengths_to_leaves);
					paralog_tree.for_each(analyse_tree);
					analyse_tree.results(sum_sqrs, avg_percent_stddev);

					args.VERBOSE() << "\tOptimize UPGMA tree using kitsch...\n";
					t_progress_indicator dots(args.VERBOSE(), 50000);
					kitsch_tree::optimize_bcpp_tree(paralog_tree, sub_cluster.names,
														sub_cluster.matrix,
														args.cnt_iterative_fitch_passes,
														args.power,
														args.keyhole_size,
														best_iteration,
														sum_sqrs,
														avg_percent_stddev,
														cnt_trees_examined,
														dots);
				}

				//888888888888888888888888888888888888888888888888888888888888888888888888

				// kitsch only

				//888888888888888888888888888888888888888888888888888888888888888888888888
				else if (args.analysis_type == eKITSCH)
				{
					args.VERBOSE() << "\tBuilding kitsch tree...\n";
					t_progress_indicator dots(args.VERBOSE(), 50000);
					kitsch_tree::kitsch_make_tree(paralog_tree,
													sub_cluster.names, sub_cluster.matrix,
													args.cnt_iterative_fitch_passes,
													args.power,
													best_iteration,
													sum_sqrs, avg_percent_stddev,
													cnt_trees_examined,
													dots);
				}


				// calculate branch lengths
				bcpp_tree_ops::do_calculate_branch_lengths branch_lengths(paralog_tree);

				// calculate sum squares
				bcpp_tree_ops::do_sum_sqrs2 analyse_tree(sub_cluster.matrix,
														branch_lengths.lengths_to_leaves);
				paralog_tree.for_each(analyse_tree);
				analyse_tree.results(sum_sqrs, avg_percent_stddev);



				//
				// Output timing and results
				//
				strms.ofs_timing << boost::format(	"%|1$14|\t"
													"%|2$15|\t"
													"%|3$11|\t"
													"%|4$14|\t"
													"%|5$9|\n")
													% sum_sqrs
													% avg_percent_stddev
													% cnt_trees_examined
													% best_iteration
													% timer.reset();

				if (!strms.printout_analysis(paralog_tree, ">Cluster# " + sub_cluster_name,
														branch_lengths,
														args.taxon1, args.taxon2))
					return 0;
			}
		}
	}




	//
	//	print clusters with missing entries
	//
	if (problems.size())
	{
		cerr << "Some matrixes had missing or duplicate pairs:\n";
		for (unsigned i = 0; i < problems.size(); ++i)
		{

			cout << format("\tCluster #%|1$4|: "
								"%|2$4| genes, "
								"%|3$4| gene pairs expected, "
								"%|4$4| gene pairs missing, "
								"%|5$4| redundant duplicates.\n")
					% problems[i].cluster_name
					% problems[i].cnt_genes
					% problems[i].cnt_expected_pairs
					% problems[i].cnt_missing_pairs
					% problems[i].cnt_duplicate_pairs;
			cout.flush();
		}
	}
	cerr << "\nAll finished!!\n\n";

	return 0;
}

